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We analyze in detail the heating of bosonic atoms in an optical lattice due to incoherent scattering 
of light from the lasers forming the lattice. Because atoms scattered into higher bands do not 
thermalize on the timescale of typical experiments, this process cannot be described by the total 
energy increase in the system alone (which is determined by single-particle effects). The heating 
instead involves an important interplay between the atomic physics of the heating process and 
the many-body physics of the state. We characterize the effects on many-body states for various 
system parameters, where we observe important differences in the heating for strongly and weakly 
interacting regimes, as well as a strong dependence on the sign of the laser detuning from the excited 
atomic state. We compute heating rates and changes to characteristic correlation functions based 
both on perturbation theory calculations, and a time-dependent calculation of the dissipative many- 
body dynamics. The latter is made possible for ID systems by combining time-dependent density 
matrix renormalization group (t-DMRG) methods with quantum trajectory techniques. 
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I. INTRODUCTION 

Recent experimental advances with ultracold quantum 
gases in optical lattices [TH3] have opened opportuni- 
ties to explore novel phenomena in many-body lattice 
physics [31-115] , including aspects of quantum magnetism 
with bosonic and fermionic atoms, and possibilities for 
characterizing the phase diagram of the Fermi-Hubbard 
model. However, the production of strongly interacting 
many-body states at the very low temperatures required 
to reach many of these phases remains a key challenge in 
current experiments [16H20] . 

In this context it is very important to be able to char- 
acterize and control heating processes arising in exper- 
iments. These can appear, e.g., via laser fluctuations 
that give rise to phase and amplitude noise on the lat- 
tice potential, through collisional losses of atoms, or via 
incoherent scattering of the lattice light. While heating 
of single atoms in dipole traps due to incoherent scatter- 
ing was characterized a long time ago [2TJ [22 , and dis- 
cussed recently in the specialized case of an optical lattice 
potential [53], heating of many-body states in strongly 
interacting systems presents a new problem due to the 
interplay between essentially single- or few-body heat- 
ing processes, and the characteristics of the many-body 
state. As a result of this interplay, different many-body 
states can be more or less sensitive to particular heating 
processes. 

Here we analyze in detail the heating of bosonic atoms 
in an optical lattice due to incoherent scattering of light 
from the lasers forming the lattice. We choose this case 
both because it provides a clear example in which to 
study the interplay between the atomic physics of the 
heating process and the many-body physics of the state, 
and because incoherent scattering is expected to be the 



dominant heating mechanism in recently analyzed exper- 
iments [16] . The Bose Hubbard model corresponding to 
bosonic atoms in the lowest-energy Bloch band of an op- 
tical lattice is defined via the Hamiltonian (ft = 1) 

H BH = -jJ24b 3 + luJ2b\ 2 h 2 , (1) 

where the first term is a kinetic energy describing the 
hopping of bosons on the lattice with amplitude J with 
hi (b\) bosonic destruction (creation) operators, and the 
second term is an onsite interaction with strength U. The 
phase diagram of the Bose Hubbard model contains a su- 
perfluid phase for J ^> U, and a Mott insulator phase for 
J <C U, connected by a quantum phase transition. This 
phase diagram has been discussed extensively in the liter- 
ature using both mean field and exact Quantum Monte 
Carlo techniques (see Ref. [16 , and references cited). 
For reference below, we note that the critical point for 
the quantum phase transition appears at (U/J) c « 3.37 
in ID [23], whereas in 2D and 3D the quantum phase 
transition at unit filling is closer to the mean field value 
of (U/(zJ)) c ~ 5.8 with z the number of nearest neigh- 
bors for each lattice site. Our goal below is to character- 
ize the effects of incoherent scattering of lattice light on 
many-body states for various system parameters. The 
non-equilibrium dynamics of this process are described 
by a master equation, which includes the coherent atomic 
dynamics in form of a multiband Bose Hubbard model, 
as well as an incoherent part describing the effects of 
spontaneous emission. In this context we observe impor- 
tant differences in the heating for strongly and weakly 
interacting regimes, as well as a strong dependence on 
the sign of the laser detuning. 

A key feature of the physics here is that it is not suf- 
ficient to determine the total rate of energy increase in 
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the system in order to characterize the heating. This is 
because some single-particle excitations do not thermal- 
ize on typical experimental timescales (e.g., individual 
atoms excited to higher Bloch bands). Indeed, the mean 
rate of energy increase is independent of the interactions 
and of the sign of the laser detuning, as is found for a 
single atom (see, e.g., Ref. [23] )■ Instead, the change in 
the many-body state must be characterized in terms of 
characteristic correlation functions for the state, e.g., the 
single particle density matrix, which characterizes off- 
diagonal order in the superfluid regime. 

We compute the time dependence of characteristic 
correlation functions based both on perturbation the- 
ory calculations, and a time-dependent calculation of 
the dissipative many-body dynamics. The time depen- 
dent dynamics are described by a many-body master 
equation, and can be computed in the mean-field limit 
using a density-matrix Gutzwiller approach. For ID 
systems they can be computed exactly by combining 
time-dependent density matrix renormalization group (t- 
DMRG) methods with quantum trajectory techniques 
[25] • We show that in the weakly interacting regime, 
bosons are strongly susceptible to heating in the sense 
that long-range order in the superfluid ground state is 
destroyed by a localization mechanism in spontaneous 
emission events. In contrast, a Mott Insulator ground 
state, in which each atoms is already exponentially lo- 
calized at a particular lattice site, is very robust against 
spontaneous emissions. The rate of destruction of long- 
range order depends on the total scattering rate, not on 
the energy input into the system, and so is much more 
rapid for red-detuned lattices than for blue-detuned lat- 
tices. 

This article is organized as follows: in Sec. [Il] we briefly 
review the description of a single atom in an optical lat- 
tice including spontaneous emission and then present the 
model we use to describe the situation of many atoms. 
In Sec. |III| we present key quantities characterizing spon- 
taneous emission such as the scattering rate, the total 
increase in energy and the effect on the key correla- 
tion functions for a system in the ground state of a the 
Bose Hubbard model. In Sec. |III 5] and |HTC] we present 
the results of fully time dependent calculations based 
on an exact numerical calculation for ID lattices, and 
a Gutzwiller mean field approach for 3D lattices, respec- 
tively. In Sec. |IV| we then present a summary and out- 
look. 



II. MODEL 

We consider bosonic atoms in an optical lattice, which 
is generated by a far-detuned laser fields via the AC-Stark 
shift, and study how spontaneous emission affects the 
many-body state of the atoms. We will describe atomic 
dynamics in terms of master equations, 



with p the reduced density operator of the atoms, where 
we trace over the bath of vacuum radiation modes, H is a 
multiband Hubbard Hamiltonian for bosons in an optical 
lattice, and £ is a Liouvillian representing the effects of 
spontaneous emission. 

In Sec. |II A] we will first write out this equation for two- 
level systems with the off-resonant excited state elimi- 
nated. Such master equations have been developed in 
the context of laser cooling [23 [23 , mainly in a single 
particle context [22, 23, 28 , and we will adapt and gen- 
eralize them to the present problem. We will also discuss 
the differences in the scattering rates and similarities in 
the heating rates for blue and red detuned laser light in 
this case. We note that for very large detunings the as- 
sumptions of a two-level atom and rotating wave approx- 
imation (RWA) break down, but we stick to this model 
because the two-level results are compact and transpar- 
ent, and are readily generalized to include more levels, 
and contributions from counterrotating terms. Finally, in 



Sec. II B we generalize to N atoms, including a discussion 
of interactions and non-idealities in typical experimental 
setups. 



A. Single particle case 

In this subsection, we briefly review the dynamics of 
a single two-level atom with mass m and internal states 
\g) and |e) in an optical field (for a more detailed dis- 
cussion of these dynamics see Ref. [22] ) . Below we begin 
from the optical Bloch equations including the atomic 
motion, and then derive an effective master equation for 
the ground state \g) in the limit of large laser detuning. 
We then derive a form for the master equation expanded 
in terms of Wannier functions for the lattice potential, 
before discussing the key features of heating in a deep 
lattice, and how this heating depends on the sign of the 
laser detuning. 



1. Optical Bloch Equations with motion 

The motion of a two-level atom driven by a laser field 
and undergoing spontaneous emission is described by Op- 
tical Bloch Equations (h = 1) [2"2"I |2"5] . 



P 



= -i[H,p}+ (2) 
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-i [H, p] + Cp, 



Here p is the reduced density matrix of the two-level atom 
{!#)> l e )}- I n Eq. j2j r is the decay rate for the excited 
state e. The Lindblad operator C u — \g){e\e~ lk " gU ' x de- 
scribes the return of the atomic electron |e) — > \g) to 
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the ground state after a photon emission in direction u, 
including the associated recoil kick to the atom. Here 
k eg = ui eg /c « fci is the wavenumber associated with the 
atomic transition frequency uj eg , and N(u) is the distri- 
bution of directions for the emitted photons. If we denote 
the unit vector along the dipole moment of the transition 
with d, this distribution is given by 



N(u) = 



8n 



1 



u ■ d 



(4) 



The atomic Hamiltonian ^ contains the kinetic energy 
and the laser interactions. The laser driving the atom has 
optical frequency ujl = ckj, with detuning A — ljl — oj eg 
from atomic resonance. The speed of light is denoted by 
c. The laser interaction is characterized by a spatially 
dependent Rabi frequency f2(x) proportional to the elec- 
tric field of the laser and the atomic dipole moment d. 
Eqs. ((2j) and ^ are written in a frame rotating with 
the laser frequency, i.e. the optical frequencies have been 
eliminated. 



2. Elimination of the excited state 

In the limit of small saturation and large detuning 
|A| 3> fi,r we can eliminate the excited state adiabati- 
cally to obtain a master equation for the external degrees 
of freedom only. Denoting the density operator for the 
motion again with p this reads: 



dt 1 



= -i(H eBP - P Hl ff ) + Jp, 



(5) 



where the non-hermitian effective Hamiltonian is given 
by 



Heff = 



P_ 

2m 

Pi 

2m 
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4A 
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Here we have identified the spatially dependent AC Stark 
shift with the optical potential V ov t{x) and ^(x) with the 
rate of light scattering. The recycling term in Eq. ^ is 
given by 



Jp = T d 2 uN{u) 
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where the operators c u {x) 



(7) 

■ x n{x)/(2A) corre- 
spond to absorption of a laser photon followed by the 
scattering of a spontaneous photon in the direction u. 
We will write the many-body master equation in a simi- 
lar form in Sec. Ill Bl 



3. Expansion in Wannier modes 

For a periodic optical potential we can expand the mas- 
ter equation in a basis of real Wannier functions that are 



exponentially localized at each lattice site (as is done 
in the standard derivation of the Bose-Hubbard model 
PQ). Here we denote the Wannier function at lattice site 
i in the Bloch band n as w^(x). Anticipating the N 
bosons case below, we introduce second quantized mode 
operators b\ that annihilate a particle at a site i in 
the band n, and obey the usual bosonic commutation 
relations. 

In the tight binding approximation, valid for a suf- 
ficiently deep lattice, we then obtain for the effective 
Hamiltonian 



n,{i,j) 

i 

~ 2 



3 Z-^ * * 
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(8) 



The first line in Eq. |8]) is a multiband Hubbard Hamilto- 
nian, and the second line are decay terms. The relevant 
matrix elements are: 



tO) 



r (n) 



d 3 xwf\x)f^ + V opt (x)\ W ( j n \x) 



I 



•r.c W ( r\x)[^-+v opt (x)) W r>(x) 



,(n), 



d 3 x (x)j(x)w\ rn ^ (x). 



(n) 

In practice the hopping rates ■ can easily calculated 
from the band structure. For the Lindblad operators in 
the recycling term we obtain 



(9) 

(10) 
(11) 



E 



(12) 

They describe the redistribution of the atoms in the 
Bloch bands due to absorption of a laser photon from 
the optical lattice lasers followed by emission of a pho- 
ton. We note that because the Wannier functions are 
exponentially localized, spontaneous emission processes 
coupling atoms between neighboring sites are small, and 
these matrix elements have been neglected above. 

In Sec. |II B| we will generalize these equations to in- 
clude inter-particle interactions. 



4- Heating processes for red and blue detuned light 

We now consider the differences in the description of 
the heating processes for red-detuned (A < 0) and blue- 
detuned (A > 0) light. These cases are distinguished by 
the sign of V op t(x), which results for the blue-detuned 
case in the minima of the potential occurring at the min- 
ima of the field intensity, i.e., the minima of fi(a;), but 
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for the minima to occur at the maximum field intensity 
in the red-detuned case. 

In this subsection we will write the potential along 
ID for simplicity, but these concepts are readily gen- 
eralized to the 3D case. In order to make a com- 
parison where a minimum of the potential is al- 
ways centered at x = 0, we write the Rabi fre- 
quency to behave as fl(x) = Qq cos(fcix) for a red- 
detuned lattice, and Cl(x) = r2osin(fcix) for a blue- 
detuned lattice. The optical potential then behaves 
as V opt (x) = (|O | 2 /l 4A l)sin 2 (fc L a:) for the blue- 
detuned case, and V opt {x) = — (|f2o| 2 /|4A|) cos 2 (kLx) — 
(\n \ 2 /\4A\)sm 2 (k L x)- (|n„| 2 /|4A|) in the red-detuned 
case. In both cases, the minima of the potential oc- 
cur at x = 0, ±Ti/k L , ±2iT/k L ■ ■ ■ = 0, ±a, ±2a .... The 
different form of for red- and blue-detuned lat- 

tices has a large effect on the type of heating pro- 
cesses that are possible, as we can see by considering 
the behavior of the field around the minima of the field 
where the atoms are trapped, where for red-detuned light 
(O re( j(x) « fi (l — (fc_Lx) 2 /2)), and for blue-detuned light 
(f2blue(a;) ~ tt k L x). 

To obtain a simple picture for the scattering pro- 
cesses, we can consider atoms that are tightly trapped 
in the lattice, so that at each site and for each dimen- 
sion, the Lamb Dicke parameter, r/ = k^a®, which com- 
pares the extension ao of the lowest band Wannier func- 
tion to the wavelength of scattered photons (2tt /k^) is a 
small parameter. The dependence of this parameter on 
the depth of the lattice V = |fi | 2 /(4|A|), is given by 
V = k L a = (W/Er)- 1 / 4 (where E R = k 2 L /(2m) is the 
recoil energy). In the limit of deep lattices (V 3> 2Sr), the 
Wannier functions can be approximated with harmonic 
oscillator wave functions around the potential minima 



and the integrals in Eq. (11) and (12) may be evaluated 
in a Lamb Dicke limit (77 <gc 1), expanding the Rabi fre- 
quency and the plane wave in a power series in kj,x. The 
leading order scattering processes in this limit are de- 
picted in Fig. [T] We then obtain that the scattering rate 
for particles initially trapped in the lowest band to lowest 
order in the Lamb Dicke parameter is given by 



r scat t = 7 (0 ' 0) 



r|Qo| 2 

4A 2 

r|Qo| 2 

4A 2 



A < 0; 
2 , A>0. 



(13) 



As is expected from the fact that atoms in the blue- 
detuned case are at the minimum of the light intensity, 
the scattering rate is reduced substantially in this case, 
by a factor of rj 2 . This corresponds to a substantial de- 
crease in the probability of atoms being scattered back to 
the lowest Bloch band in this case, where in the Harmonic 
oscillator approximation the rate of scattering returning 
atoms to the lowest band is given by 



pj.b. , 

scatt 



r|a | 2 

4A 2 „ 

r|n | 



A < 0; 



4A 2 'I 1 *-* u - 



der term for the blue-detuned lattice is scattering into the 
first excited band, which is of order rj 2 . In the Harmonic 
oscillator approximation this coupling to the first excited 
band occurs at the identical rate in the red-detuned lat- 
tice . 

While the red-detuned and blue-detuned cases exhibit 
considerably different scattering rates, the rate at which 
the energy of the atom E = (H) increases is identical 
for red and blue detuning, and completely independent 
of the initial motional state of the atom (see Appendix 
|A), 



■ _ r|o | 2 
e -^a 2 ~ Er - 



(14) 



Whereas scattering into the lowest band is the dominant 
term (order rf) in the red-detuned lattice, the leading or- 



In the harmonic oscillator approximation, the indepen- 
dence of this rate on the sign of the detuning can be seen 
clearly, because an increase in energy arises solely from 
coupling to higher bands, and whilst the rate of scat- 
tering events returning particles to the lowest band are 
very different for blue-detuned and red-detuned light, the 
rate with which particles are scattered to higher bands is 
identical in the two cases. When the same computation 
is performed with Wannier functions, the lowest Bloch 
band has a finite energy width, and thus a small energy 
increase is obtained from atoms scattered into the low- 
est band. The heating rate in the red- and blue- detuned 
cases remains identical, however, with the scattering rate 
into higher bands in the red-detuned case being slightly 
reduced to compensate the energy increase from scatter- 
ing into the lowest band. 

We note at this point that the rate of energy increase 
in the system does not completely characterize the heat- 
ing process, and we will show below that for many-body 
systems, heating in blue- and red- detuned lattices can 
function very differently, despite the fact that the rate of 
energy increase remains independent of the sign of the 
detuning. A key physical characteristic of spontaneous 
emission processes when characterizing how the many- 
body state changes is that they tend to localize particles 
undergoing the scattering event on a lengthscale of the 
wavelength of the emitted photon. This will be described 
m more detail in Sec. Ill CI 



5. Remarks on the master equation 

We conclude with a few remarks regarding the validity 
of the above model. The assumption of a two-level sys- 
tem here is clearly an oversimplification in some respects. 
However, it is very clear how to generalize these results 
to more realistic atomic models and experimental setups. 

First, in order to generate an isotropic 3D cubic opti- 
cal lattice, it is typical to use three laser beams that are 
independent, either because they are slightly detuned or 
because they have orthogonal polarizations. The opti- 
cal potential is a sum of optical potentials for all three 
beams, and all terms in the master equation must be 
summed over contributions from the different beams. In 
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FIG. 1: Schematic picture of the leading-order heating pro- 
cesses via scattering of laser photons in a red-detuned (left) 
and a blue-detuned (right) lattice, in the Lamb-Dicke limit 
(see text). In a red-detuned lattice the fastest process returns 
the scattered atom to the lowest band, whilst this process is 
strongly suppressed in the blue detuned lattice. In both cases 
there are higher order processes in the Lamb-Dicke parameter 
r\ of the lattice, in which the atoms are scattered into higher 
Bloch bands. In the limit of a deep lattice, these processes 
occur at identical rates in the red-detuned and blue-detuned 



the calculations below we explicitly choose different lat- 
tice depths in different directions, and perform this sum- 
mation, in order to produce results that are as close as 
possible to current experiments. The sum of the effec- 
tive scattering rate |f2o| 2 r/(4A 2 ) over the beams will be 
denoted 70. 

Second, for typical parameters, the lattice lasers are 
so far detuned that the rotating wave approximation is 
not strictly valid. However, in the limit where the excited 
state can be eliminated, this does not change the effective 
description of the physics except for small quantitative 
modifications in the prefactor of the scattering rate and 
the optical potential. 

Finally, in multi-level atoms, the far-detuned lasers will 
couple to many excited states. Again, in the limit of large 
detuning, all of these states can be adiabatically elimi- 
nated to produce the same effective model with prefactors 
that arise from correctly summing the contributions from 
all excited levels. 

Similar remarks also apply to the iV-Atom Master 
equation that we present below. 



B. TV- Atom Case 

We now present the full model for N atoms in the lat- 
tice. Below we first state the full master equation assum- 
ing two-level atoms, and discuss the origins of each term, 
including terms arising from interactions and collisional 
loss in the system. 



1. Master equation for N atoms 

For N bosonic atoms the evolution of the reduced 
system density operator p is given by a master equa- 
tion which we write in using second quantization. Elim- 
inating the atoms in the excited state and defining 
field operators ip(x) with bosonic commutation relations 

i){x), ^(y) = S(x — y) for atoms in the ground state, 

we derive again a master equation of the form 



H cS p 



pHts 



Jp, 



with non-hermitian effective Hamiltonian 

rad 1 zjxoll 



tr tt 1 rrrad , zircon 



(15) 



(16) 



The first contribution to the effective Hamiltonian, Hq, is 
the term describing motion of single atoms in the optical 
lattice, 



H = J d 3 xft(x) (- 



2m 



+ V opt (x) i/>(x), 



(17) 



which is the same form derived in Sec. Ill Al 

The radiative part of the master equation describing 
the couplings of the atoms to the vacuum modes of the 
electromagnetic field are contained in the effective Hamil- 
tonian 



//;,;r = / / d 3 xd 3 y m{y ^* {x) G(k eg (x - v))ft(*)1>Hv)${v)fo) 
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4A 2 



F{k eg {x - y))^ \x)^ '{y)i>{y)i>{x) 



(18) 
(19) 



and recycling term 



F(k eg (x - y))ft (x)tp(x)p^ (y)i>(l 



(20) 
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with functions F and G denned as 
F(i) = J d 2 uN(u)e- lu -t 

= I{^H<M)°) + (l-3(<M?)(^-^)}, (21) 
1 f°° dC C z 

^{-( 1 -(^.fl^ + (l-3(rf- £ Y)(^ + ^)}. (22) 



where V denotes the principal value integral. 

Finally, we have an effective Hamiltonian accounting 
for short range collision physics in the presence of laser 
fields, as well as associated losses 



rrcoll 



il '.v ( g(x) — i^7 2 (a;) j $ (x)ffi (x)ip(x)iip(x) 



(23) 

where the functions g and 72 will be discussed below. 

We will now give a more detailed discussion of the ra- 
diative and collisional contributions to the master equa- 
tion, as well as commenting on the validity of these equa- 
tions in each case. 

Radiative terms - As mentioned above, radiative pro- 
cesses are contained in the effective Hamiltonian H^ d 
and the recycling term J . These contributions to the 
master equation are obtained by eliminating the vacuum 
modes of the radiation field for an ensemble of two-level 
atoms, as discussed by Lehmberg [501131], followed by an 
adiabatic elimination of the excited state. The present 
master equation is a straightforward generalization of 
Lehmberg's .ZV-atom Bloch equations by including the 
quantized motion of the atoms, and by writing the master 
equation in second quantized form. The first line in H^ d 
is the dipole-dipole interaction due to exchange of pho- 
tons between the atoms. The second line contains a single 
particle decay term corresponding to the absorption of a 
laser photon, followed by spontaneous emission, as dis- 



cussed in Sec. II A The second term in Eq. (19) is a col 



lective radiative term associated with super- and subra- 
diance. The functions F and G appearing in Eqs. ( 18|19 ) 
and (20) are defined in Eqs. ( 21|22 1. For distances much 



smaller than the optical wavelength £ = fc e9 |x — y| <C 1, 
i.e. for particles on the same lattice site, the function G 
approaches the static dipole-dipole interaction diverging 
as r -3 , while F(0) = 1. For distances larger than the 
wavelength, i.e. particles on distant lattice sites, both G 
and F fall off in an oscillatory manner on a lengthscale 
set by the wavelength of the emitted photons. A plot of 
can be found in Fig. [2] 
We conclude the discussion of the radiative terms with 



three remarks. First, the recycling term J in Eq. (20) in 



Thus a spontaneous emission event will localize a parti- 
cle within a wavelength. This is a key mechanism behind 
destruction of long-range order due to spontaneous emis- 
sion processes, and we will return to this discussion below 



volves Lindblad operators in the form of atomic densities 
^ (x)ij)(x) smeared out by the function F(k eg (x — y)). 



in Sec. |IIC| and Sec. |III A 4| Second, it is easily checked 
that the effective Hamiltonian together with the 

recycling term ( 20 ) give a trace preserving master equa- 
tion in the usual Lindblad form. This follows from the 
commutation relations of 4>{x) together with F(Q) = 1. 
Thirdly, we note that the dipolar function G(£) must 
be understood as regularized on a short distance scale 
of molecular interactions where the two-level model un- 
derlying the derivation of the radiative master equation 
breaks down. In typical experimental situations, where 
the lattice is very far detuned, the dipole-dipole inter- 
action term will give only very small corrections to the 
effects of normal collisional interactions, which we dis- 
cuss in the next paragraph. In the derivation of a Hub- 
bard model, these dipole-dipole interaction terms will 
give small corrections to the onsite interaction energy, 
and will be absorbed into this value for the purposes of 
the discussions in following sections. 

Collisional terms — We now turn to the collisional 
terms describing short range scattering. In the absence of 
laser light, and for scattering at sufficiently low energies 
and densities, this short-range physics provides a bound- 
ary condition for the two-body scattering wavefunctions 
at longer distances, so that the resulting dynamics can 
be described as an effective two-body interaction with a 
single parameter, the scattering length a s [32Tl34j . Thus, 
this short range collision physics is accounted for by a 
contact potential in the man y p article Hamiltonian, as 
represented by if^ 11 in Eq. (23) with g — <iTrh 2 a s /m. 
The scattering length a s can be calculated by solving a 
set of coupled channel equations for scattering of atoms 
|35j . It is, at least on a conceptual level, straightforward 
to include in these channels not only collisional processes 
but also light assisted collisional interactions. An ex- 
ample is provided by the discussion of the optical Fes- 
hbach resonances by Fedichev et al. [31)H3"8] . where for 
red detuned laser light excited electronic states of the 
atom-atom complex provide resonances and thus a reso- 
nant enhancement of the ground state scattering length. 
Such a light-modified scattering length will reflect the 
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local laser intensity, and thus, in principle, be spatially 
dependent. We account for this by writing an intensity- 
dependent and spatially varying contact coupling g(x) 
in Eq. ( 23 1 . Away from the optical Feshbach resonance 
we expect <?(x) ~ 4irh 2 a s /m. More important, radiative 
loss in these processes can be accounted for by an imagi- 
nary part of the a scattering length, which will again be 
intensity and thus spatially dependent, as discussed in 
the context of optical Feshbach resonances. We account 
for this loss by including an intensity dependent imag- 
inary contribution to the scattering length —172 in the 
contact interaction ( p3| . Along a similar reasoning losses 
in collisions between chemically reactive molecules have 
recently been modeled by imaginary scattering length. 

In general, there can be corrections to collisional in- 
teractions due to higher order processes, e.g., three-body 
collisions that can generate three-body losses. These are 



neglected here. Also, in the discussions below we would 
like to separate the role of spontaneous emission events 
from that of light-assisted collisions. As a result, we will 
set 72 = for the purposes of calculations exploring the 
effects of heating due to spontaneous emissions. 



2. The master equation in a Wannier basis 

Multi-band Hubbard model - When the potential cor- 
responds to an optical lattice, we obtain a multi-band 
Bose-Hubbard model from Eq. ( 15 ) for the coherent part 



of the evolution by expanding the field operators in a 

Wannier basis, ip(x) — Eni 1 "! W^l"'' an< ^ a PPly m g 
the assumption of local tunneling and interactions in a 
deep lattice pfl, e.g. for an isotropic 3D lattice, 



(24) 



Aid) 



with tunneling rates j\ n ^ a nd onsite interaction energy shifts J7( fe > J > m > n ) arising from the contact and dipole-dipole 



interactions in Eqs. (18 



23 1. Below we will denote j\ n ^ and JJ^ n,n,Tl,n ^ for the lowest band as J and U, respectively. 



The remaining part 01 the master equation in this basis takes on a similar form to that in Sec. |II A| 

P = ~i [H,p] + dp. 

The term describing scattering of laser photons, denoted dPi can be written in this basis as 



dp 



1 



-5V 



(k,l,m,n) 



b (fc)t 6 (0 



Am) ti(n) 



2 

Here the sum runs over fc, I, m, n} and the matrix elements for the different processes are defined as: 
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(fe,j,- 



_ 



d 3 xd 3 y 



4A2 



F{k eg (x - y))wf\x)wf\x)w^ n '(y)wf'(y). 



(25) 



(26) 



(27) 



Single-band Hubbard model - For red-detuning, where 
the dominant processes in the scattering will return 
atoms to the lowest band, it can be convenient to focus 
on the physics resulting from these dominant processes. 
If, in addition, we neglect the small terms that arc not di- 
agonal in a position basis (~ F(k e , g a)), then ( 15 ) reduces 
to 

p = -i[H B H, p}+^2 lijiipui - l/2nimp - l/2pn l n i ). 

(28) 

Here 7 is the effective scattering rate, and Hbh is the 
Bose-Hubbard Hamiltonian for the lowest band, as given 
in Eq. |l| with b, = ^ 0) , n, = b\b h and e = = 0. 
This form of the master equation is used below for quan- 
tum trajectories simulations, and also clearly demon- 
strates the localizing effects of spontaneous emissions 



that are discussed in the following subsection. The ex- 
pressions here for U and J are also used when describing 
the initial state of the system, where atoms are assumed 
to be confined to the lowest Bloch band. 



Localization of atoms due to spontaneous 
emission events 



As mentioned above, a key characteristic of sponta- 
neous emission events, both for single particle and many- 
particle systems, is that they localize the atom undergo- 
ing the scattering event on a length scale given by the 
wavelength of the emitted photon 39J . This is most 
clearly seen from the recycling term in Eq. ( 20 ) , where 
the function F(k eg (x — y)) determines the coherence be- 
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FIG. 2: Plot of the function F(£) for £ parallel to d (see 
Eq. (211 for the definition of F). Note that the localized 



form of this function reflections the localization of atoms in 
spontaneous emission events. 



tween emissions at different points in space, x and y. 
In Fig. [2] we plot the function F(£), showing that it is 
clearly localized. The length scale of the localization is 
then determined by l/k eg , which is proportional to the 
wavelength of emitted light. When we write the master 
equation in a Wannier basis, the fact that F(k eg (x — y)) 
is localized on a length scale similar to the lattice period 
means that the master equation will be approximately 
diagonal in position space, with the recycling term given 
by 



(k,l,m,n), (fe) f r (J) i(m)ti(n) 
Ti.i b i b i P b i b i ■ 



This diagonal form in b\bi will dephase coherent super- 
positions in which atoms are delocalized over many sites 
into mixtures in which for each classical possibility the 
atom is localized. This is a representation of the fact that 
information about the atom's position is transmitted to 
the environment by the emitted photon, and in the ab- 
sence of a measurement of the photon, we are left with a 
classical mixture of different possible locations where the 
photon could have been scattered. Below we will see that 
this mechanism is very important when one considers the 
effect of spontaneous emissions on many-body states. In 
particular, this localization tends to destroy long-range 
order, which is a key property of superfluid states. 



many body states in perturbation theory as described by 
the master equation (15). The calculations are based on 



the full multi-band master equation, but use perturba- 
tion theory in the limit of small 7/ J, j/U, where 7 is the 
effective scattering rate, and also neglect interactions of 
atoms that have undergone spontaneous emission events. 
We quantify heating for different many-body states, and 
for different signs of the laser detuning. We first investi- 
gate the effective scattering rate, and the rate of energy 
increase, and find results that are very similar to the 
single-particle case presented in Sec. |II A[ but which are 
modified quantitatively due to collective effects. How- 
ever, because atoms scattered to higher bands will typ- 
ically not thermalize with those in the lowest band on 
experimental timescales, the increase in mean energy is 
not sufficient to determine the change in the many-body 
state due to heating. Instead, it is important to directly 
investigate changes in the characteristic correlation func- 
tions. We investigate the change in correlation functions 
in subsections MI A 41 and IIIIB1 



In Sec. IIIB we combine quantum trajectories tech- 
niques with t-DMRG methods in order to account for 
longer time evolutions and thermalization after sponta- 
neous emission events by propagating the master equa- 
tion directly in time. However, for these purposes we find 
it convenient to restrict calculations to the lowest band 



master equation (28), relevant for deep red detuned lat- 
tices. This gives us an exact solution, but only for one 
dimensional lattices. 



Finally in subsection III C we present a mean field de- 
scription based on Gutzwillcr mean field theory to in- 
vestigate the master equation (15) for a 3D lattice. We 



expect this approach to provide a semi-quantitative de- 
scription for the heating processes, but this will certainly 
not capture details for the spatial dependence of the first 
order correlation functions, and their time dependence 
due to heating. 



A. Perturbation Theory 



III. QUANTIFYING DECOHERENCE AND 
HEATING 



1. Scattering Rate 



We start our discussion of decoherence and heating in 
subsection |III A| by studying the evolution of different 



We can determine the scattering rate for spontaneously 
emitted photons directly from Eq. (15), and we find 
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r 
r 



l xd?yF{ke 9 {x ~ y)Mx)n(y)Tr {ft (x)^(x)pft (y)^(y)} 



d A xQ(x 



) 2 Tr{ft(x)ftx)p} 



r 

4A2 



d 3 xd 3 yF(k eg (x - y))n(x)fl(y)Ti {ft (y)ft (x)fty)ftx)p} 



(29) 



As in the single particle case discussed in Sec. |II A 4[ this 
rate depends strongly on the sign of the detuning, and 
it now also depends strongly on the characteristic many- 
body state through the correlation functions computed 
for p. In Fig. [3]we plot the scattering rate for atoms con- 
fined to move in ID by an anisotropic lattice, for which we 
can compute the many-body ground state of the Hubbard 
model using t-DMRG methods, and hence the correlation 
functions that appear in Eq. (29). The scattering rate is 



plotted as a function of U/J for a system at unit filling, 
and we note immediately that the scattering rate is again 
much larger for red-detuned light A < than for blue- 
detuned light A > 0. The scattering the red-detuned 
case is again dominated by processes returning atoms to 
the lowest Bloch band, as shown by the dashed lines in 
Fig. [3] In the superfluid regime (i.e. for U/J < 3.37), 
the rate for processes that leave the atom in the low- 
est band is increased by bosonic enhancement |40) . due 
to the presence of other delocalized atoms at the place 
where the spontaneous emission event occurs. Since the 
process that returns an atom to the lowest band is much 
more prominent for red-detuned light, this enhancement 
is much more visible in that case than for blue-detuned 
light. As the interaction strength is increased, atoms 
tend to be exponentially localized on individual sites in 
the Mott-Insulator phase, corresponding to U/J> 3.37, 
and for increasing U/J this enhancement is strongly re- 
duced. 



1.5 



0.5 



A< 



A> 



10 
U/J 



15 



20 



FIG. 3: (color online) Comparison between the total process 
rate r sca tt/7o (solid lines) and the rate for processes back to 
the lowest band (dashed lines) for A < (upper lines) and a 
A > (lower lines) in the ground state of a ID lattice (Using 
DMRG ground states for Bosons in a ID lattice: V x = V v = 
30 Er, V z = 10 -Eh) The rates include scattering from all three 
lattice generating beams. 



many-body state. Furthermore, the increase in energy per 
particle is the same as in the single particle case presented 

TTKa 



Sec. 



However, as we will discuss in the next 
subsection, this is not the key quantity for determining 
changes to the many-body state. 



2. Total rate of energy increase 



3. Thermalization and atoms in excited Bloch bands 



From the master equation Eq. (15), we can calculate 



the rate of change d(H) /dt of the total mean energy of 
the atom cloud due to scattering of laser photons. Here 
H is the hermitian part of the effective Hamiltonian in 
equation ( 15 ), and C\p is the part of the master equations 



that describes incoherent scattering of the laser photons: 



~{H)=Tr{H£ lP }. 



For each standing waves with laser wavenumber hi , giving 
an effective Rabi frequency Q,(x) = cos(fcjx) we find 
(see Appendix |A| 



dt 



(H) 



rjOoj 

4A 2 



E R N. 



(31) 



Thus, the total rate of energy increase is independent 
of the sign of the detuning, and of the properties of the 



The dominant contribution to the energy computed in 
Sec. Ill A 2 is due to atoms scattered into higher Bloch 
bands, which can easily be shown in the Lamb Dicke 



limit (see Sec. II A 4). The energy gain in a scattering 
event in which the atom remains in the lowest band is on 
the order of the tunneling rate J, whereas atoms scat- 
tered into higher bands gain an energy greater than the 



(30) bandgap energy £ gap . However, because e gap 3> J, U, it is 



not possible for these atoms that have been scattered to 
higher bands to thermalize this energy input on typical 
experimental timescales. Indeed, it would require very 
high-order processes in perturbation theory in J/e gap and 
U/e g ap in order to return these atoms to the lowest band. 
As a result, these processes that create a large change in 
the mean energy of the state will nonetheless often con- 
tribute very little to a change in the many-body state. 
Thus it is clear that the total energy does not completely 
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FIG. 4: (Color online) Relative rate of change of the 
(integrated) single particle density matrix 5115(21,22) = 
ff dxdy(ffl {x,y,zi)^){x,y,Z2)) in an effectively ID lattice 
(14 = V y = 30 E R ; V z = WEr). In the tightly bound 
transversal directions the atoms are assumed to be in the 
lowest band. Scattering from all three lattice generating laser 
beams is taken into account (with weights corresponding to 
the lattice depths). (The lattice constant is denoted by a.) 



characterize the change in the many-body state. Instead, 
we must look directly at other quantities, such as corre- 
lation functions determining the characteristic properties 
of the many-body state, and how they are changed as a 
result of the heating processes. We will show below that 
although the rate of energy increase is independent of the 
sign of the detuning, that the change in character of the 
many-body state can be strongly dependent on that sign. 



Single particle density matrix and momentum 
distribution 



For bosons in an optical lattice, a key quantity in char- 
acterizing the many body state is the single particle den- 
sity matrix S(x 1 ,x 2 ) = (£Ci)i/j(a? 2 )) • In the super- 
fluid regime, the system exhibits off-diagonal long range 
order (or quasi- long range order for ID systems) , whereas 



in the Mott Insulator regime this function decays expo- 
nentially as a function of distance \x± — x%\. In first or- 
der perturbation theory in 70/ J, Jo/U, the change of the 
single particle density matrix due to scattering of laser 
photons is given by 



S(x!,x 2 ) =Tr{ft( Xl )4>{x 2 )£ lP } 

1 r 

~~ 24A2 
-2F(k(x! 



(32) 

n( Xl ) 2 + n(x 2 ) 2 
x 2 ))n( Xl )n(x 2 )){ft(xS(x2)}. 

(33) 



Here we have also neglected interactions with atoms that 
have undergone spontaneous emissions. These will be 
treated in Sec. IIIB Noting that F(k(xi — x 2 )) — > for 
\xi — a?2 1 a the off-diagonal long range order changes 
according to: 



\x 1 —x 2 \^>a 



1 r 

24A 5 



(n( Xl ) 2 + n(x 2 ) 2 ) (ft{ Xl )i)(x 2 )). (34) 



The relative rates of change for the single particle den- 
sity matrix as a function of x^ and x 2 are plotted in 
Fig. [4] for red and blue detuned light, normalized to 70. 
We see a clear difference between the results for red and 
blue detuning. The reason for this is that the breakdown 
of long-range correlations is rooted in the localization ef- 
fect of spontaneous emission events, as was discussed in 
Sec. |II C[ This localization depends not on the rate of 
energy input into the system, but rather on the total 
scattering rate, which is larger in the red-detuned case. 
Hence, the breakdown occurs substantially faster for the 
red-detuned lattice than for blue detuning. 

At the same time, the localization effect is much more 
harmful for a superfluid state in the weakly interacting 
limit, than for a Mott Insulator state, where the particles 
are already exponentially localized at different sites. In 
the extreme limit U/J — > 00, the only significant change 
in the state is the (relatively rare) transfer of some atoms 
to higher bands. 

Note that the single particle density matrix is also 
directly connected to the momentum distribution n(p), 
meaning that the measurements that are made in exper- 
iments (including the comparison made with quantum 
monte-carlo calculations in Ref. [16] ) are directly related 
to the changes in these correlation functions. Specifically, 



n(p) = 



1 



(2tt) 3 



d 3 x 1 d 3 x 2 S{ Xll x 2 )e~ lp{xi - X2) , (35) 



so that the presence of (quasi-) off-diagonal long range 
order in the superfluid regime gives rise to peaks at re- 
ciprocal lattice vectors. The breakdown of long-range 
order due to spontaneous emission events is then directly 
related to a decrease in the visibility of these peaks, as 
we show in Fig. [5]). Again, in the extreme Mott Insu- 
lator limit, the state is very insensitive to spontaneous 
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FIG. 5: (Color online) Momentum distribution (in arbitrary 
units) of the ground state of noninteracting particles in a ID 
lattice of depth V = WEr (black), and after a time 7ot = 0.5 
if the lattice is generated by a red laser (red, dashed). The 
inset shows the central peak initially (black, solid) and after 
7o£ = 0.5 in a red (red, dashed) and in a blue lattice (blue, 
dashed-doted). (36 particles on 36 sites) 



emissions, with the initial momentum distribution be- 
ing given by the Fourier transform of the lowest band 
Wannier function, and changing only due to processes 
where an atom is scattered into higher bands. Since the 
timescales for these processes are similar for red and blue 
detuning, both lattices have similar effect on the Mott 
Insulator. 

Using perturbation theory calculations it is difficult to 
quantify the effects of interactions between particles on 
the lattice after the spontaneous emission has occurred, 
in particular effects due to partial thermalization of the 
energy added to the system. In order to address this, it 
is necessary to find a method to propagate the master 
equation directly in time. This is discussed in the next 
section. 



B. Propagation of the single band master equation 
§28^ with Quantum Trajectories & TEBD 



In order to further quantify the heating process, espe- 
cially the effects of interactions in the system after spon- 
taneous emission events have occurred, we now inves- 
tigate means to directly integrate the master equation 
(15 1. For ID systems we can compute time-dependent 



expectation values from the master equation exactly by 
combining time-dependent density matrix renormaliza- 
tion group (t-DMRG) methods with quantum trajectory 
techniques, as discussed in Ref. [25 . In order to sim- 
plify the numerical computation, we focus on processes 
returning atoms to the lowest Bloch band, which are 
dominant for red detuning, as discussed above, and base 
the simulations on the single-band master equation given 
in Eq. (28). In the 2D/3D case, it is numerically pro- 



hibitively expensive to apply these methods directly for 
realistic system sizes, however we expect that certain as- 
pects can be described semi-quantitatively by a mean- 
field theory treatment similar to that in Ref. [36], which 
will be discussed in Sec. MI CI 

In the combination of quantum trajectories and t- 
DMRG, t-DMRG [3TrH5] provides a convenient means 
to propagate states that are not too far from equi- 
librium in a ID system, whereas quantum trajectories 
[44l [45] is a method to compute time-dependent cor- 
relation functions from the master equation based on 
a stochastic propagation of states. Briefly, the idea is 
that each stochastic trajectory begins from an initial 
pure state (sampled from the initial density matrix) , and 
is propagated based on the non-Hermitian Hamiltonian 
Hbh ~ (7/2) TiiUi, except for at randomly 



QT 
eft 



sampled times iy, where quantum jumps occur, 



(36) 



corresponding to the localization of a particle onsite due 
to the spontaneous emission event. In the stochastic 
simulation the times tj are points where the norm of 
the state falls below a randomly chosen threshold. At 
these times, a random site L is selected for the spon- 
taneous emission event according to the probabilities 
Pi j oc ( , 4'(tj)\ n i- n ij\' l l J ('tj)) an d applied. Expectation val- 
ues are computed based on stochastic averages, which 
converge rapidly as a function of the number of trajecto- 
ries. Each of the computations presented here was based 
on sampling 1000 trajectories, which was possible on a 
timescale of a few days with a cluster of ca. 100 CPU 
cores. In every case, numerical convergence was checked, 
with the number of states retained for decompositions in 
t-DMRG being computed to x = 200. 

Results from these calculations are shown in Figs. [6j 
[7] and [§} In Fig. [6] we show the increase in energy 
for a system with 36 particles on 36 sites interaction 
strength U/J — 3 and 7 = 0.01J. We observe ex- 
cellent agreement with results from perturbation the- 
ory at short times, where we have the linear relation 
(H BH ) t fa (H BH ) + f t (H BH )\ t=0 t, with the rate of 
change of the mean energy from Eq. _|28|) given by 

we similarly 
peak in the 



(H 



iJEuMb,)- In Fig. 7 



see the characteristic decay of the centra, 
momentum distribution as off-diagonal correlations de- 
cay over time. 

In Fig. H we show the evolution of the off-diagonal 
elements of the correlation functions S(i,j) = (b\bj), 
where in contrast to the results presented in the previ- 
ous section, we now include the effects of collisions after 
spontaneous emission events, i.e., partial thermalization 
of the energy added to the system. The basic change to 
these previous results is immediately apparent in Fig. [8^,, 
which shows off-diagonal correlations in the superfluid 
regime for UJ J = 2 and in the Mott Insulator regime for 
U I J = 10 as a function of time. Where the perturbation 
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FIG. 6: Development of the mean energy in a system with 
36 particles on 36 sites, at U/J — 3 just below the super- 
fluid to Mott Insulator transition and 7 = 0.01J, computed 
from Eq. (28 1 via quantum trajectories methods, compared 



with the equivalent result from first order perturbation the- 
ory. Statistical error bars are shown based on the stochastic 
average over 1000 trajectories. 




FIG. 7: Momentum distribution (in arbitrary units) calcu- 
lated for a ID lattice (36 particles on 36 sites) with U/J = 3 
(i.e. in the superfluid regime) and a scattering rate 7 = 0.01J 
from Eq. ( 28 1 via quantum trajectories methods in a lattice of 
depth V — IOEr. Averages are taken over 1000 trajectories. 
The scale of the momenta axis is 1/a, where a is the lattice 
constant. 



theory results, neglecting scattering after a decay event, 
give us identical rates of relative decay of these values 
for different interaction strenghts, here the rates depend 
strongly on U/J. To further quantify this, we plot values 
of the correlation function normalized to the same initial 
value in Fig. [8}j, for different values of U / J ranging be- 
tween U/J — 2 and U/J = 10. From perturbation the- 
ory calculations not including collisions between atoms 
after spontaneous emission events, these should all de- 
cay at the same rate, as ji{b\bj) — —j(b\bj). However, 
we see that in the superfluid case, the rates decay more 
rapidly at intermediate times, as thermalization creates 
further decay of the off-diagonal correlations in addition 
to the initial localization effect of the spontaneous emis- 
sion events. In contrast, for the Mott-Insulator limit, the 
results deviate already strongly from the perturbation 
theory results on a timescale given by 1/t/, and on longer 
timescales the off-diagonal elements are barely changed, 




FIG. 8: Comparison of the decay of off-diagonal correlations 
in the single particle density matrix S(i,j) = (b\bj) as a func- 
tion of time. These are computed from the master equation 
for the lowest band Eq. ([28]) with j/J = 0.01 for 24 parti- 
cles on a 24 site lattice by combining quantum trajectories 
methods with t-DMRG. (a) Comparison of the decay of off- 
diagonal elements £\ S(i,i + x) when U/J = 2 (superfluid 
regime, upper lines) and U/J — 10 (Mott Insulator regime, 
lower lines). In each case we show correlations at separation 
distances x = 1 (solid lines) and x = 2 (dashed lines), each 
of which are averaged over i. We see that the correlations for 
the Mott Insulator state at U/J = 10 are almost constant, 
whereas for U/J = 2, the decay becomes more rapid than at 
initial times, (b) Comparison of S(i,i + 2) averaged over i 
and normalized to the value at time t = 0, \S(i, i + 2)[[t = 0] 
for U/J = 2,4,6,8, 10 (seen here from bottom to top). The 
dashed line shows the corresponding result from perturbation 
theory. Each computation was averaged over 1000 trajecto- 
ries, and error bars are shown in (b). For (a), the statistical 
errors fit inside the line thickness. 



as the strong interactions tend to result in the small local 
correlations being reestablished. It is worthwhile to note 
that longer-range correlations in the Mott-Insulator state 
can actually slightly increase in comparison with their 
initially exponentially small values, as energy added to 
the system is thermalized. 
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Propagation of the master equation via a 
mean-field ansatz 



Guzwiller ansatz assumes a product state over the lattice 
sites, 



In two or three dimensions, exact computation of time- 
dependent dynamics becomes numerically extremely ex- 
pensive for all but very small system sizes. However, it is 
possible to gain insight from the time-evolution based on 
a mean- field ansatz. Below we employ a treatment of the 
master equation analogous to the Gutzwiller mean-field 
analysis of the ground states of interacting systems (valid 
in higher dimensions), where we generalize these ideas to 
time-dependent dynamics. Such a Gutzwiller mean field 
approach has been developed in Ref. to describe dy- 
namical quantum phase transitions as a competition be- 
tween coherent Hamiltonian and incoherent Liouvillian 
dynamics. This method is simple to apply, and allows us 
to easily include higher bands in the computation, which 
can be numerically expensive for the exact methods dis- 
cussed in the previous section. A simple picture of the 
heating can be built up in terms of the distribution of par- 
ticles in different bands and the increase in energy and 
entropy in the system, as well as decay of the superfluid 
order parameter due to heating, if we begin in a super- 
fluid state. The inclusion of higher bands here also allows 
us to quantitatively see the lack of thermalization of en- 
ergy input as atoms transferred to higher bands. As with 
ground state calculations based on a typical Gutzwiller 
ansatz, the results are expected to be more accurate for 
systems in higher dimensions, and it is not possible to ex- 
tract accurate information about the spatial dependence 
of correlation functions, which we have already discussed 
in some detail for the ID case. 

For zero temperature, ground state calculations of 
the single band Bose Hubbard Hamiltonian Eq. ([I]), the 



|*)=n^) = nE/. 



Wl 



(37) 



The amplitudes /„ of states with occupation number n 
at site i are used as variational parameters to minimize 
(Hbh) with the constraints of normalized wave function 
and a fixed mean number of particles per lattice site. 
Here we identify the mean field Hamiltonian, 



H 



MF 



E 



j 



(V4 



i. 



+ m) + -Ub} b: 



as a sum of local operators. The superfluid phase 
manifests itself in a nontrivial superposition |</>j) = 
with a nonzero expectation value of the de- 



h\bj\ 



whereas in the 



J2n fn 

struction operator ipi = ^2 
Mott Insulator phase with filling uq one has p n = S n ^ no , 
such that ipi = 0. At unit filling this transition from su- 
perfluid to Mott phase occurs at U/(zJ) rj 5.8, where z 
denotes the number of nearest neighbors. 

In the same spirit we use a factorization ansatz to solve 
our master equation ( 24 )— ( 27 1: we write p = <^> i Pi with 
Pi = Tr^\{p} [46], and derive an equation of motion for 
the density operator at site i: pi = Tr^;{£p}. Then the 
multiband equation of motion is given by 

pi = -i [flMF,i, Pi] + £-MF,iPi- (38) 

The multiband mean-field Hamiltonian at lattice site i 



huF,i= E M? OfM n) + (bf)b^) + e w&w>) + e \u^™^bf^bfh^ ] b { r\ (39) 



n.m.k.l 



where the brackets denote the expectation value under the density operator p: (...) = Tr{. . . p}. The incoherent part 
at lattice sites i is: 



C 



n - ST Ak,l,m,n) / ,(fc)tr (I) n i(m)fL(n) 
MF,iPi ~ 2^ 7 M I "i °i Pi "i "i 



n,m,k,l 



i&wt^/oyo l - Pi b {m) h {n) b {k) h {l) 



(40) 



In particular, for a homogeneous situation pi = a, the 
above master equation is a nonlinear equation for a, since 
the Hamiltonian depends on the expectation values of the 
destruction operators in the various bands. We note that 
in contrast to the problem studied in }46) the dissipative 
part is linear in p. This comes from neglecting scattering 
processes in which the atoms change lattice site and and 



, (k,l,m,n) 

the symmetry 7^ ^ 



7. 



■ifc.l) 



It is easy to show 



that the master equation ( 38 1 preserves the trace and the 
mean particle number. 

Here we will consider the example of homogeneous sys- 
tem (pi = a for all sites) in an isotropic cubic 3D lattice 
(z = 6), and will restrict our discussion to the lowest 
band and first three excited bands (which are degenerate 
for the isotropic case). In this way we include the dom- 
inant heating channels in our calculations, as in a rela- 
tively deep lattice the Lamb Dicke parameter 77 is small, 
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FIG. 9: Particle number statistics for the lowest Band and one 
of the three first excited bands, computed for a noninteracting 
system as a function of time using a Gutzwiller ansatz for the 
system density operator. Initially the statistics in the lowest 
band are Poissonian, corresponding to a coherent state. The 
statistics in the lowest band change towards an exponential 
decay with no as a function of time, while the population 
in the exited band increased. (Calculated with 7 (°' ' > > = 
O.OlzJ in a isotropic 3D lattice (z — 6) of depth V = WEr.) 
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FIG. 11: Time evolution of the mean number of particles 
per lattice site in the lowest and first excited bands, com- 
puted using a Gutzwiller ansatz for the system density opera- 
tor, and plotted as a function of time for different interaction 
strengths. As the rate of population of higher bands is inde- 
pendent of the interaction strength, we see from the similarity 
of results for interacting and non- interacting systems that pro- 
cesses returning the atoms to the lowest band are very slow, 
leading to a lack of thermalization of energy transferred to 
the system by transferring particles to higher bands. As for 
Fig. [9] we choose 7 (°' ' '°> = O.OlzJ in a isotropic 3D (z = 6) 
lattice of depth V = 10E R . 




FIG. 10: Particle number statistics for the lowest Band and 
one of the three first excited bands, computed for an inter- 
acting system with U/(zJ) = 1 as a function of time using 
a Gutzwiller ansatz for the system density operator. As in 
Fig. [9] we choose 7 (°' ' >°) = O.OlzJ in a isotropic 3D lattice 
(z = 6) of depth V = WE R . 



and the rates of heating to higher bands are significantly 
smaller, as discussed in Sec. |II A| We begin each calcu- 
lation with the Gutzwiller ground state at zero temper- 
ature for atoms in the lowest band of the lattice, with a 
given U I J value at unit filling. The state is then propa- 
gated in time using the above method, with a Gutzwiller 
ansatz for the system density operator. The scattering 
rates are based on a red-detuned optical lattice with a 
depth of V = 10E R , setting 7 (°'°> ' ) = O.OlzJ. 

In Fig. [9] and Fig. [10] we show the evolution of the 
particle number distribution in the lowest and excited 
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FIG. 12: Evolution of the square modulus of the superfluid or- 
der parameter |i/>(t)| 2 , computed using a Gutzwiller ansatz for 
the system density operator under the Eq. ( |38[ ) for different 
interaction parameters. The initial states are the correspond- 
ing Gutzwiller ground states. The different initial states have 
different non- vanishing superfluid order parameters, which in 
all cases decays on the timescale set by the localization rate 
in the lowest band. (Calculated with 7 (°> ' '°> = O.OlzJ in a 
isotropic 3D lattice (z = 6) of depth V = WEr.) 



bands as a function of time for a non-interacting gas 
(U = 0), and an interacting gas with U/{zJ) = 1. For 
vanishing interaction (U/(zJ) — 0) the initial state ex- 
hibits Poissonian number statistics in the lowest band, 
P(n ) = (n |o-|n ) = (n } no e-<™°> /(n \). We start with 
one particle per lattice site in the lowest band so that 
(no) = 1- As the state evolves, we see that the num- 
ber statistics change rapidly from Poissonian towards an 
exponential distribution. At the same time, the prob- 
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FIG. 13: Evolution of the mean total energy per lattice site, 
computed using a Gutzwiller ansatz for the system density op- 
erator. The initial values depend on the interaction strength 
of the system. The increase in energy is the same in all cases, 
as we saw in Sec. |III A2| Since we have included only one 
excited band, the total energy saturates, leading to the de- 
viation from the linear increase of the exact solution. The 
thin dashed lines in dicates the increase in energy as calcu- 
lated in Sec. 



Ill A 2 



(Calculated with 7 (°' ' ' ) = Q.QlzJ in 
a isotropic 3D lattice [z = 6) of depth V = 10-E.r.) 



co 



2.5 
2 

1.5 
1 

0.5 




/✓ y 

/ / s 
fi s 
ft / 
ft / 




// / 
/// 


U/(zJ)=0 


h 


- - - U/(zJ)=1 


p 


U/(zJ)=3 







100 



200 
zJt 



300 



400 



FIG. 14: Increase of entropy per lattice site for different inter- 
actions, computed using a Gutzwiller ansatz for the system 
density operator. Starting in the ground state, the entropy 
is initially zero. Due to the different initial states the in- 
crease in entropy is higher for lower interaction parameters, 
reflecting the more significant change in more weakly inter- 
acting states due to spontaneous emission events. (Calculated 
with 7 (°' ' '°) = O.OlzJ in a isotropic 3D lattice of depth 
V = WE R .) 



In Fig. 11 we show the time evolution of the mean par- 
ticle number for the lowest and excited bands for different 
interaction strengths. We note that for different interac- 
tion strengths, these numbers are identical up to very 
long times. As the rate of population of higher bands is 
independent of the interaction strength, we see from the 
similarity of these results that that processes returning 
the atoms to the lowest band via collisions are very slow, 
leading to a lack of complete thermalization of the energy 
transferred to the system, as was discussed in Sec. |III A~3] 
In Figs. [l2"}{T4| we quantify the heating for states of dif- 
ferent initial U/(zJ) values via different quantities. The 
quantity corresponding to off diagonal long range order 
is = ">P an d we identify this quantity as an 

order parameter for superfluidity In Fig. 12 we show 
the evolution of this superfiuid order parameter, where 
we clearly see the destruction of superfluidity for atoms 
in the lowest band as a result of the heating processes. 
This is analogous to the destruction of long-range order 
discussed in |III A4[ In the Mott Insulator phase, the su- 

and 

In Fig " 



pcrfluid order parameter is zero at the beginning 
remains zero throughout the evolution. 
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we 

show the total energy increase in the system over time. 
These results agree well for short times with the results 
obtained in Sec. IIII A 21 A nice feature of the Gutzwiller 
ansatz is that it is simple also to calculate the increase 
in entropy for the system, as we start from a pure state 
and heat the system into a mixed state. The entropy per 
lattice site is plotted in Fig. fT4l and describes the same 



basic behavior as Figs. 12 and [13) 

In summary, we find that a product ansatz for the sys- 
tem density operator in the spirit of a Gutzwiller mean- 
field treatment gives a simple semi-quantitative picture 
for the heating process. Including higher Bloch bands, 
we observe quantitatively that particles heated to higher 
bands are not transferred back to the lowest band on 
typical experimental timescales, even in the presence of 
significant interactions. As in the previous section, we see 
that scattering in the lowest band gives rise to a destruc- 
tion of superfluidity in the system, here characterized via 
the superfiuid order parameter. This is always zero in 
the Mott Insulator phase, and the key properties of this 
phase in the Gutzwiller description change only in that 
particles can be heated to higher bands. In this sense, we 
see that more strongly interacting states are significantly 
more robust against heating due to spontaneous emission 
events. 



ability of population in the higher bands gradually in- 
creases. Similar behavior is observed when the inter- 
action strength is increased to U/(zJ) = 1. When we 
choose interaction strengths above the superfluid-Mott 
Insulator transition at unit filling, which is not shown 
here, the dynamics of the number statistics becomes triv- 
ial. In the Gutzwiller representation, we have exactly one 
particle per lattice site, and the only dynamics arising 
from spontaneous emission events is gradual population 
of excited bands. 



IV. SUMMARY AND OUTLOOK 

We have shown that the heating of atoms in optical 
lattices due to spontaneous emission events is strongly 
dependent on the characteristics of the lattice, especially 
the detuning of the lattice beams, and also on the many- 
body characteristics of the state. Because atoms scat- 
tered to higher Bloch bands will not thermalize with 
atoms remaining in the lowest band on experimental 
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timescales, it is not sufficient to compute the total rate 
of energy increase in order to determine the change in 
the many-body state. Instead, the heating can be char- 
acterized, e.g., by computing characteristic correlation 
functions, such as the single particle density matrix for 
bosons. 

We found that the higher scattering rate for red- 
detuned lattices as opposed to blue-detuned optical lat- 
tices corresponds to a much more rapid breakdown in 
off-diagonal order in a superfluid state, due to the local- 
izing effect of spontaneous emission events. In contrast, a 
Mott Insulator state, where the atoms are already expo- 
nentially localized, can be strongly robust against spon- 
taneous emission events. 

In an experiment, other design considerations will have 
to be taken into account when choosing the lattice de- 
tuning, e.g., the different rates of light-assisted collisions 
72 for red- and blue-detuned laser light. This interplay 
is particularly interesting because light-assisted collisions 
tend to be more prominent for blue-detuned light, which 
is where the rate of spontaneous emissions is lowest. For 
production of states where the atoms are exponentially 
localized at different lattice sites, red detuned lasers could 
be used without strong adverse effects. However, for pro- 
duction of states with off-diagonal long-range order, we 
have shown here that the laser detuning is an important 
consideration. In the future, the quantum trajectories 
methods we have used here could be extended to include 
two-body loss terms, and used to analyze the competi- 
tion between different heating mechanisms in these ex- 
periments. 

Another key future direction will be the investigation 
of heating of fermionic species. The results here give 
an indication that states in which atoms are localized, 
e.g., a Mott Insulator state with possible additional spin- 
ordering could, under favorable circumstances, be rela- 
tively robust against spontaneous emission events. 



DARPA OLE program. 
Appendix A: Calculation of the total heating rate 



The total heating rate (31 1 is calculated as the change 



rate of H , the hermitian part of H c g in Eq. ( 15 1 
d 



dt 



[H) = Tr{H£ lP } = 



i r 

24A7 



Tr{ J J d 3 xd 3 yF(k(x - y))fl(x)fl(y) 
H,ft(x)Tp(x)],ft(yMy)]p}. (Al) 



We use the approximation k eg ss k^ = k and introduce 
the notation T(x,y) = F(k(x — j/))fi(x)0(y). Noting 
that F(0) = 1;VF| = and AF| = -fc 2 and using 
Maxwell's equation A£l(x) = — k 2 fl(x) we readily derive 
the following relations: 

A x T(x, y)\ y=x = -2k 2 n(x) 2 , (A2) 
V x T(x, y)\ v=x = Q{x)V x Q{x). (A3) 

Here the A denotes the Laplacian. Further we have: 



[H,ft(x)$(x)] 
-1 
2m 



A, 



^(x)^(x)-^(x)(A x i>(x)^y (A4) 



Using this together with the relations (A2) and (A3) in 



Eq. (All we find after partial integration: 
d 1 



(H) 



Tr{ / d 3 x (2k 2 fl(x) 2 ft (x)ip(x)+ 



dt x "' 2 2mA 2 

-(v x n(x) • v x n(x) + n(x)A x n{x))i> 1i (x)^(x)J p}. 

(A5) 
(A6) 
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In the lattice we have the relation 2k 2 fl(x) 2 + W x ft(x) 
V x Cl(x) + Q(x)A x Cl(x) = k 2 , which leads to: 



ml k 

4A 2 2~r 



-N. 



(A7) 



This is the result stated in Eq. (31 ) 
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